
rm(list = ls())

###################################################
###################################################
###################################################

# Load Packages

pacman::p_load('readxl', 'dplyr', 'data.table', 'ggplot2', 'ggthemes', 'ggview', 'patchwork', 'showtext', 
               'fixest', 'broom', 'kableExtra', 'sf', 'fastDummies', 'rdrobust', 'tibble')

###################################################
###################################################
###################################################

# Baseline Panel

df <- read_excel("data/reg_panel_final.xlsx") %>% data.table()

df[, ln_purge := log(purge+1)]
df[, popu57 := popu[year == 1957], by = countyid]
df[, ln_purge_pc := log((purge+1)/popu57*1000)]
df[, over_exc := ifelse(retain < 165, 1, 0)]
df[, ln_dr := log(dr)]
df[, glf := ifelse(year %in% 1958:1961, 1, 0)]

df[, dr_1957 := dr[year == 1957], by = countyid]
df[, popu57 := popu57 / 10^6, by = countyid]
df[, ln_clan := log(genealogy+1)]
df[, ln_ccp := log(ccp_pc_intp)]
df[, ln_kmt := log(kmt_pc_intp)]

###################################################
###################################################
###################################################

# Baseline Covariates

psc <- c('FM', 'AM', 'age', 'party', 'edu', 'longmarch')
loc <- which(names(df) %in% psc)
names(df)[loc] <- paste0('psc_', names(df)[loc])

cov_cnty <- c('popu57', 'weather', 'grain_pc57', 'school', 'ln_clan', 'num_cc', 'ln_ccp', 'ln_kmt')
cov <- c(cov_cnty, names(df)[loc])

###################################################
###################################################
###################################################

# RDD Data

rd <- read_excel("data/RD_data_final.xlsx", guess_max = 5000) %>% data.table()

rd[, grid := paste(floor(long), floor(lati), sep = "-")] # cluster by 1-degree grids

rdfe <- cbind(dummy_cols(floor(rd$long), remove_first_dummy = T, remove_selected_columns = T), 
              dummy_cols(floor(rd$lati), remove_first_dummy = T, remove_selected_columns = T),
              dummy_cols(rd$EPROV, remove_first_dummy = T, remove_selected_columns = T)) # longitude, latitude, and province FEs

rd[, dist := ifelse(south == 1, dist_abs, (-1)*dist_abs)]
rd[, ln_purge := log(purge+1)]
rd[, ln_dr := log(dr_1960)]
rd[, ln_popu57 := log(popu57)]
rd[, ln_clan := log(genealogy+1)]
rd[, ln_ccp := log(ccp_pc_intp)]
rd[, ln_kmt := log(kmt_pc_intp)]
rd[, suit_rice := suit_rice / 100]
rd[, rug_index := rug_index / 100]
rd[, rain := rain / 100]

v_a <- c('ln_ccp', 'ln_kmt', 'num_cc')
v_b <- c('dr_1957', 'ln_popu57', 'grain_pc57', 'school', 'ln_clan', 'dist_rail')
v_c <- c('suit_rice', 'rug_index', 'rain', 'temp')
rd_var <- c(v_a, v_b, v_c)

###################################################
###################################################
###################################################

# Variable Labels

lab <- c('sat_all' = 'All satellites',
         'sat_grn' = 'Grain satellites',
         'over_exc' = 'Over-extraction',
         'ln_dr' = 'Log mortality rate',
         'log(cohort)' = 'Log birth cohort sizes',
         'glf' = 'GLF',
         'purge' = 'Purge victims',
         'purge_pc' = 'Victim density',
         'ln_purge' = 'Log purge victims',
         'ln_purge_pc' = 'Log victim density',
         'log(purge_intp+1)' = 'Log purge victims (interpolated)',
         'popu57' = 'Population, 1957',
         'ln_popu57' = 'Log population, 1957',
         'weather' = 'Rainfall deviation', 
         'grain_pc57' = 'Grain productivity', 
         'school' = 'Middle schools', 
         'ln_clan' = 'Log genealogies', 
         'num_cc' = 'Central elites', 
         'ln_ccp' = 'Log CCP density', 
         'ln_kmt' = 'Log KMT density',
         'dr_1957' = 'Death rate, 1957', 
         'cps_native_p' = 'Native party secretary',
         'cps_ht_north' = 'Northern party secretary',
         'log(zf_num_kill+1)' = 'Log suppression victims',
         'log(anticor)' = 'Log anticorruption victims',
         'countyid' = 'County',
         'year' = 'Year',
         'dist_rail' = 'Distance to railroads',
         'suit_wheat' = 'Wheat suitability',
         'suit_rice' = 'Rice suitability',
         'rug_index' = 'Terrain ruggedness',
         'rain' = 'Precipitation',
         'temp' = 'Temperature')

###################################################
###################################################
###################################################

# Figure Layout

mytheme <- function(){
  theme_minimal() +
  theme(axis.title.x = element_text(margin = margin(t = 10)),
        axis.title.y = element_text(margin = margin(r = 10)),
        axis.text = element_text(size = 11),
        strip.text = element_text(size = 12),
        plot.title = element_text(hjust = 0.5),
        strip.background = element_rect(fill = 'gray85', colour = NA))
}

update_geom_defaults("hline", list(linewidth = .3, linetype = 2))
update_geom_defaults("vline", list(linewidth = .3, linetype = 2))

###################################################
###################################################
###################################################

# Regression Table Layout

setFixest_etable(style.tex = style.tex(main = "aer", yesNo = c("Yes", "No"), fixef.suffix = " FEs"), 
                 digits = 'r3', digits.stats = 2, fitstat = c("n", "ar2"))

###################################################
###################################################
###################################################


